Example: Blowfly Model

Likelihood-free inference for the blow-fly model was introduced by Simon N. Wood. We model here the discrete time stochastic dynamics of the size $N$ of an adult blowfly population as given in section 1.2.3 of the supplementary information.

where $e_t$ and $\epsilon_t$ are independent Gamma random deviates with mean 1 and variance $\sigma_p^2$ and $\sigma_d^2$, respectively.

using Distributions, StatsBase, LikelihoodfreeInference
Base.@kwdef struct BlowFlyModel
    burnin::Int = 50
    T::Int = 1000
end
function (m::BlowFlyModel)(P, N₀, σd, σp, τ, δ)
    p1 = Gamma(1/σp^2, σp^2)
    p2 = Gamma(1/σd^2, σd^2)
    T = m.T + m.burnin + τ
    N = fill(180., T)
    for t in τ+1:T-1
        N[t+1] = P * N[t-τ] * exp(-N[t-τ]/N₀)*rand(p1) + N[t]*exp(-δ*rand(p2))
    end
    N[end-m.T+1:end]
end

Let us plot four realizations from this model with the same parameters.

using StatsPlots
plotly()
m = BlowFlyModel()
plot([plot(m(29, 260, .6, .3, 7, .2),
           xlabel = "t", ylabel = "N", legend = false) for _ in 1:4]...,
     layout = (2, 2))

To compare different realizations we will use histogram summary statistics. In the literature one finds also other summary statistics for this data.

summary_statistics(N) = fit(Histogram, N, 140:16:16140).weights

summary_statistics (generic function with 1 method)

We will use a normal prior on log-transformed parameters.

function parameter(logparams)
    lP, lN₀, lσd, lσp, lτ, lδ = logparams
    (P = round(exp(2 + 2lP)),
    N₀ = round(exp(4 + .5lN₀)),
    σd = exp(-.5 + lσd),
    σp = exp(-.5 + lσp),
    τ = round(Int, max(1, min(500, exp(2 + lτ)))),
    δ = exp(-1 + .4lδ))
end
(m::BlowFlyModel)(logparams) = m(parameter(logparams)...)
target(m::BlowFlyModel) = [(log(29) - 2)/2,
                           (log(260) - 4)*2,
                           log(.6) + .5,
                           log(.3) + .5,
                           log(7) - 2,
                           (log(.2) + 1)/.4]
lower(m::BlowFlyModel) = fill(-5., 6)
upper(m::BlowFlyModel) = fill(5., 6)
prior = TruncatedMultivariateNormal(zeros(6), ones(6),
                                    lower = lower(m), upper = upper(m))

TruncatedMultivariateNormal{Distributions.MvNormal{Float64,PDMats.PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},Float64}(
mvnormal: DiagNormal(
dim: 6
μ: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
Σ: [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0]
)

lower: [-5.0, -5.0, -5.0, -5.0, -5.0, -5.0]
upper: [5.0, 5.0, 5.0, 5.0, 5.0, 5.0]
)

Let us now generate some target data.

model = BlowFlyModel()
x0 = target(model)
data = summary_statistics(model(x0))

1000-element Array{Int64,1}:
 0
 0
 0
 0
 1
 0
 0
 5
 1
 1
 ⋮
 0
 0
 0
 0
 0
 0
 0
 0
 0

Adaptive SMC

smc = AdaptiveSMC(prior = prior)
result = run!(smc, x -> summary_statistics(model(x)), data,
              maxfevals = 2*10^5, verbose = false)
using PrettyTables
pretty_table([[keys(parameter(zeros(6)))...] quantile(smc, .05) median(smc) mean(smc) x0 quantile(smc, .95)],
             ["names", "5%", "median", "mean", "actual", "95%"],
             formatter = ft_printf("%10.3f"))

┌───────┬────────────┬────────────┬────────────┬────────────┬────────────┐
│ names │         5% │     median │       mean │     actual │        95% │
├───────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│     P │     -2.194 │      0.507 │      0.035 │      0.684 │      1.980 │
│    N₀ │     -1.709 │     -0.036 │      0.008 │      3.121 │      1.732 │
│    σd │     -1.583 │     -0.008 │     -0.060 │     -0.011 │      1.551 │
│    σp │     -1.774 │     -0.069 │      0.020 │     -0.704 │      1.990 │
│     τ │     -1.259 │      0.316 │      0.259 │     -0.054 │      1.816 │
│     δ │     -1.483 │      0.141 │      0.114 │     -1.524 │      1.747 │
└───────┴────────────┴────────────┴────────────┴────────────┴────────────┘

histogram(smc)

corrplot(smc)

KernelABC

k = KernelABC(prior = prior, delta = 1e-1, K = 10^3, kernel = Kernel())
result = run!(k, x -> summary_statistics(model(x)), data)
pretty_table([[keys(parameter(zeros(6)))...] quantile(k, .05) median(k) mean(k) x0 quantile(k, .95)],
             ["names", "5%", "median", "mean", "actual", "95%"],
             formatter = ft_printf("%10.3f"))

┌───────┬────────────┬────────────┬────────────┬────────────┬────────────┐
│ names │         5% │     median │       mean │     actual │        95% │
├───────┼────────────┼────────────┼────────────┼────────────┼────────────┤
│     P │     -1.884 │      0.426 │      0.157 │      0.684 │      1.923 │
│    N₀ │     -1.680 │      0.025 │      0.032 │      3.121 │      1.852 │
│    σd │     -1.744 │     -0.016 │      0.000 │     -0.011 │      1.727 │
│    σp │     -1.589 │     -0.030 │      0.031 │     -0.704 │      1.775 │
│     τ │     -1.548 │      0.123 │      0.079 │     -0.054 │      1.713 │
│     δ │     -1.572 │      0.046 │      0.001 │     -1.524 │      1.655 │
└───────┴────────────┴────────────┴────────────┴────────────┴────────────┘

histogram(k)

Kernel Recursive ABC (with callback)

k = KernelRecursiveABC(prior = prior,
                       K = 100,
                       delta = 1e-3,
                       kernel = Kernel(bandwidth = Bandwidth(heuristic = MedianHeuristic(2^3))),
                       kernelx = Kernel());

We will use a callback here to show how the estimated parameters evolves.

using LinearAlgebra
res_krabc = run!(k, x -> summary_statistics(model(x)), data,
                 maxfevals = 1300,
                 verbose = true,
                 callback = () -> @show norm(k.theta - x0)/norm(x0))

(x = [0.9591272318875804, 0.17411703667175185, -0.2878025875347092, -0.0011779811585626268, 0.3928382643808126, 0.1929694375869359],)

This page was generated using Literate.jl.